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Abstract 


We utilize the homotopy analysis method to find eigenvalues of fractional 
Sturm—Liouville problems. Inasmuch as very few papers have been devoted to 
estimating cigenvalues of these kind of problems, this work enjoys a particular 
significance in many different branches of science. The convergence of the 
homotopy analysis method is also considered on the high order fractional 
Sturm-—Liouville problem. The numerical results acknowledge the ability of 
the proposed method. Eigenvalues are computed within a couple of minutes 
CPU time at core i3, 2.7 GHz PC. 
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1 Introduction 


In this paper, the following special class of nth order fractional Sturm— 
Liouville eigenvalue problems are considered: 
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S- gj(2)D*y(x) + go(x)y(x) = Aw(a)y(x), asa <b, (1) 
j=l 
where w(x) and q;(x), 7 = 0,1,2,...,n, are integrable functions over [a, 0]. 


The (left-sided) Caputo fractional derivative D°’ of order a; € (j — 1,9), 
j =0,1,2,...,m, is defined by 


ee © (A(t) 
DW) = Tena f, Gone 


Note that n should be an even number denoted by 2r. The separate boundary 
conditions of problem (1) are as follows: 


y (a) = 0, (2) 
y)(b) = 0, (3) 


fori € S’ C S := {0,1,2,...,2r — 1} in which S$’ has r members. Problems 
(1), (2), and (3) generally have arisen from linear fractional equations (see 
[3, 4, 14, 8]). 


Definition 1. The left-sided Riemann—Liouville fractional integral opera- 
tor of order a is defined by 


Ila) = 55 | “(@ — 1) ly(t)at, (4) 


where y € L'[0,T] and a € Rt. 

Some useful properties of the operator J® are found in [18, 19, 16]. It must 
be mentioned that the left-sided Caputo fractional derivative (4) is originally 
defined via the left-sided Riemann-Liouville fractional integral (4) [12], as 
follows 

D%y(x) = IP (x) = y(z), > 0, 


where a € Rt, n= [a] +1 and y € L'(0,7). 

Fractional Sturm—Liouville problems (FSLP) are interesting problems 
from a practical point of view. It turns out that such fractional equations have 
seen many applications in science and engineering problems such as captur- 
ing the dynamical behaviors of amorphous materials, for example, polymer 
and porous media [15] or the superior modeling of the anomalous diffusion in 
materials with memory, for instant, viscoelastic materials for which the mean 
square variance grows faster (superdiffusion) or slower (subdiffusion) than in 
a Gaussian process, see [9]. The interested reader may wish to examine the 
great variety of works on the subject such as [1, 2, 7, 5, 6] 

Therefore, many scientists resort to mathematicians to deal with this kind 
of problem. The homotopy analysis method (HAM), Adomian decomposition 
method (see [11, 4, 14, 20]), fractional differential transform method (see [10]), 


Using homotopy analysis Mmethod to find the eigenvalues ... 51 


iterative approximation method (see [17]), and variational method (see [12]) 
were implemented on second-order FSLP. 

To the best of our knowledge, few studies have been carried out to find the 
eigenvalues of high order FSLP (see [11]). So, we generalize the homotopy 
analysis method on this particular problem. By taking polynomials as basis 
functions, we approximate the homotopy series solution that satisfies initial 
conditions (2). Then, by imposing boundary conditions (3) on this series 
solution, a system of equations will appear for which its determinant should 
be zero to have a nontrivial solution. In spite of using fewer terms to find 
the eigenvalues of the problem (1), we gain more accurate results compared 
with [11]. Moreover, our method is easily applied in solving an arbitrary 
high-order fractional Sturm—Liouville problems. On the other hand, the con- 
vergence of the series solution of the HAM for high order FSLP is proved 
which demonstrates the efficiency of this proposed method. 

This paper is organized as follows: In Section 2 some basic definitions 
of fractional calculus are presented. The HAM is mentioned, in detail, in 
Section 3. The convergence of the HAM on high order FSLP is established in 
Section 4. The numerical results are presented in Section 5 to illustrate the 
efficiency of the proposed method. The last section is devoted to including 
discussions and conclusions. 


2 Homotopy analysis method 


To illustrate the basic concept of the HAM, we consider the following general 
nonlinear functional operator: 


N[u(x)] = 0, (5) 


where N is a nonlinear operator, x denotes an independent variable, and u(x) 
is an unknown function. For the sake of simplicity, we ignore all boundary 
or initial conditions, which can be treated similarly. Through generalizing 
the traditional homotopy method Liao in [13], the following the so-called 


zero-order deformation equation was constructed: 
(1 — p)L|6(2; p) — uo(x)] = ph (x) N|6(2; p)], (6) 


where p € [0,1] is an embedding parameter, h # 0 is a nonzero auxiliary 
parameter, H(a) is an auxiliary function, Z is an auxiliary linear operator, 
uo(a) is an initial guess of u(x), and ¢(2;p) is an unknown function. It is 
important, that one has a great deal of freedom to choose auxiliary operators 
in the HAM. Obviously, when p = 0 and p = 1, the following relations, 
respectively, hold: 


$(x;0) =uo(x), (2; 1) = u(z). 
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Thus, as p increases from 0 to 1, the solution ¢(z;p) varies from the initial 
guess to the solution u(a). Expanding ¢(a;p) in Taylor series with respect 
to p, one has 


(a; p) = uo(w) + YS) um(x)p™, (7) 
m=1 
where 1 a"d(a:p) 
: 7 mr x; p 
Um (2) a ml! ay Lees 

If the auxiliary linear operator, the initial guess, the auxiliary parameter h, 
and the auxiliary function are properly chosen, series (7) converges at p = 1, 
and one has 


+00 
u(2) = ug (x) + .S Um(2), (8) 


which must be one of the solutions of the original nonlinear equation, as 
proved by Liao (see [13]). The governing equation can be deduced from the 
zero-order deformation equation. Define the vector 


Um = {uo(x),ui(2),...,Un(x)}. 


Differentiating (6), m times with respect to the embedding parameter p and 
then setting p = 0 and finally dividing them by m!, we have the so-called 
mth order deformation equation 


L[Um(2) — XmtUm—1(a)] = RH (2) Rm(U@m—1); (9) 
_ 1 a"™IN[6(e:p)] 
Rm (trait) = a ei (10) 
and 
a eT (1) 
In equation (1) we choose the initial guess as follows: 
ug (@) = ey py (@) + Copa(x) + +++ + Crpr(x) 
where j;(x),7 = 1,2,...,7r, are chosen in such a way that uo(x) satisfies the 


boundary condition (2). By implementing the HAM we can consider N first 
terms of series (8) 


N 
wn (2) = u(t) + $5 Um(z), 


as an approximation of the solution of equation (1). It is worth to note that, 
according to our main problem (1) including the parameter » and process 
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of the HAM consisting of the parameter fi, the approximate solution wy (z) 
enjoys these two parameters, that is, \ and fi. In fact , we have 


wy (a; A, h) = uo(a) + » Um (x; A, h). (12) 


By imposing the boundary condition (3) on the approximate solution (12), 
we get 


N 
wh (b; A, h) = uf? (b) + S> uO (6; d, A), 
m=1 


where i € S’. Therefore, we encounter with a system of r equations and r 
unknown coefficients ¢), C2,...,¢p. In order to have a nontrivial solution, the 
determinant of this system of equation including parameters » and A should 
be equal to zero. So, by plotting this implicit function of A and h, so-called 
h-curve, we can determine eigenvalues that correspond to each horizontal 
plateaus. 


3 Convergence of the HAM in FSLP 


In the following theorem, we consider the convergence of the HAM on frac- 
tional Sturm—Liouvill problem. 
Theorem 1. If the following series 


co 


Z(x) = uo(x) + a Um (2), (13) 


m=1 
is convergent in which u,,(x) is obtained from (9), (10), and (11), then the 


series (13) can be a solution of (1) and (2). 


Proof. Since series (13) is convergent, we have 


lim tm(x) = 0, (14) 


moo 
by the necessary condition of convergent series. 


Now, by using (11), we have 


3 [um (x) — Xmtm—1(x)] = ur(x) + (ua(x) — ui (x)) +++» + (un (x) — Un—1(2)) 


m=1 


= Un(z). 


Then by virtue of (14), we get 
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noo nm—- Co 


pa [umn (a =~ X¥mUm— 1(x)| = lim [bg (@) — Nptlenaa (a) = lim up(x) = 0. 


According to above relation and linearity of the auxiliary operator L, we can 


write 


£( ¥-fum(2) ~ mtmn-a(2))) = ye 1(tum(2) ~ xmtim-a(@)}) = 


m=1 = 


So, by means of equation (9), we have 


> 1([um() —_ Xntin-a(0)]) = AH (x) LZ [Rin (t@m—1)| = 0. 
m=1 m=1 
For as much as, h 4 0 and H(x) 4 0, the following equality is obtained 
x [Rin (Wm 1) =0 (15) 
m=1 


Consequently, from (10), one can write 


n 


Fel aed = ae [Sowier “om “| + qo(x (=)(>U uj(x)p"| 


m! Op™—- ‘pym—L 
aw(aQe ah 
= Lae [ttm—1(#)] + Go(#)tm—1 (a) — Aw(z2)tm—1 (2). 


Then, 
> Rm[@m-—1] = > [S5aj(z)D% [tm—1] + qo(x)Um—1 — Aw(x)um—1] 
m=1 m 
= Vue ood Um—1 +o. um—1] — Aw(z)[} > um—1] 
m=1 


= $0 aj(e)DZ(2)] + g0[Z(x)] — Aw(x)[Z(2)).- (16) 
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Moreover, from the initial guess uo(x) and the initial condition (2), we obtain 


u(a) =u(a) =0, ES’, 


which implies 


Z(a) = uo(a) + S$) Um(a) = uo(a) = 0. 


Therefore, Z(x) is an analytic solution of the differential equation (1). 


Remark 1. We consider {111 , 13? |. 7?7~1¥ 10m} in which r € N is 
an arbitrary constant, as a basis function of the HAM in problem (1) in which 
Qy = A2 = An_1 = 0, with subject to initial condition (2) for which even order 
of derivatives are just assumed. Then, we can write the series solution (8) 
as a linear combination of basis functions with coefficients aj;,1, @i,2,..., Qi,r, 
1=0,...,m, as follows: 


co m 
5 5 Cr mada oh Gar Ae ates Giga Pe 


m=0 i=0 
t-l ™m 
= lim y y beige ape Or a ne ge 
too 
m=0 i=0 
= lim (a9.12 + apo2? +--+ +49,-27"~1) 
t4o00. ; : 
Qr-1 1+on Qr-1tan 
+ (@o1% + +> app??? + ayaa TO fs bay pa) 
Qr-1 1+(t-Lan 
“4 (ao,12 fe eee] a0,r& . Se + Qt-1,1% +( Je 


Bk ate alles gg es) 


= lim t(ao,1.2 + apgx° + +++ + a9,-2"7"~") 
too : i > 


+ (t= 1Wlarre!te fo aye) 
14+(t-l)ay +4 Be brit Less) 


tee t (Qa 
t-1 
= lim Sct a R\lapatt fee yp age, 


too 
k=0 


ob dete 


In order to establish the convergence of above series we should have 


t—-k—-1 plt+(k+lan 
lim ( i )an+1,5% <1 
too] (t— kag, gatt (Ron) 
; ; : ; t—k-1 
for 7 =1,...,r. It is obvious that the lim; oe se <1 for all k <t. 
So, CR+LI pom <i. 


ak, j 


56 J. Biazar, M. Dehghan and T. Houlari 


4 Numerical examples 


To illustrate the proposed approach, some examples are presented in this 
section. We use Mathematica software to calculate eigenvalues of following 
problems. 


Example 1. Consider the following fourth-order FSLP: 
D*ly(x)] = Ay(2), #2 (0,1), (17) 
where 3 < a < 4, subject to the boundary conditions 
y(0) = y"(0) = 0, (18) 


and 
y(1) = y"(1) =0. (19) 
We assume that the solution of (17) can be expressed by the sct of the 


following basis functions: 


‘ 4 i 1 1 
ee a cack 


According to the boundary condition (18), we can choose the initial ap- 
proximation of the solution in the form uo(xz) = ca+ c,z°. It is obvious that 
one must choose the auxiliary linear operator 


L[$(x;p)] = D* ola; Pp). 
From (17), we define the following nonlinear operator: 
N|6(#;p)] = D° (a; p) — A0(@}P). (20) 
We have the zero order deformation equation (6) with the initial conditions 
G(0;p) =0, 9" (0;p) = 0. 
From (10) and (20), we have 


1 OD" (uo + pur +++) = Alu + pun +++) 
(m—1)! Op™-1 
= D°tUm_1(2) — AUm_-1(z). 


Riles) = 


Now, the solution of the m—th-order deformation equation (9), for m > 1 
becomes 


Um(£) = XmUm—1(£) + AJ°[H (x) Rm (te m—1)]- 


According to the rule of solution expression denoted by (8) and from (9), the 
auxiliary function H(z) is uniquely determined with H(x) = 1. 


Using homotopy analysis Mmethod to find the eigenvalues ... 57 


The first few terms of HAM series solution for a = 3.7 are as follows: 


ur = ci[(1 + h)x — 0.013788hAx +] + c2[(1 + A)ax? — 0.002166hAx>+%], 
ug = e1[(1 + h)?x + (—0.027757 — 0.0275757hAw1*+% + 0.000010h? \2a1+2%)] 
+c2[(1 + fi?)a> + (—0.004332 — 0.004332h) fix? +% + 6.40592 x 1077? 203429], 


Consequently the N—th order approximate solution of the HAM, wy (2), is 


in the form 
N 


wn (x) = Um(x) = c1¢n (2) + conn (a). (21) 


m=0 
By imposing the boundary condition (19) on (21), we get 


c1Gn (1; A, h) + conn (1; A, hf) = 
erCny(1; A, B) + camy (1; A, A) = 


In order to have nontrivial eigenfunctions, we just need to solve 


Cn (134, ) nw (13 A, h)\ 
( W(1;A,h) nll (1; d,h)) ~ ° 


The above equation, which depends on X and fh, is used to plot the h-curve. 


We can observe that in the plot of \ as a function of h, several horizontal 
plateaus occur, each of which corresponds to an eigenvalue of the problem. 


The results of Table 1 has been obtained by choosing h = —1. The exact 
eigenvalues of this problem for the integer case a = 4, are known \, = (k7)*, 
k > 1 (see [11]). The obtained results of HAM are compared with those 
calculated in [11]. As we can observe in Figure 2 the nth eigenfunction has 
exactly, n — 1 zeros in the mentioned interval. 


Table 1: Comparison results of first four eigenvalues of Example 1 


Xe a=37 a=39 @ = 3.9999 a=4 
HAM fi] HAM [i] HAM fi] HAM {] 
y 91412292 91.412293 93.533242 93.533230 97 404401 97.404401 97409091 97.409091 
de 944.796194  944.795695 | 1324.156313 1324.156357 | 1558.290698 —1558.290698 | 1558.545456 —1558.545456 
ds 4544.318681  4544.336950 | 6456.132806 6456.132485 | 7888.515088 7888.515134 | 7890.136374 — 7890.136374 
M 12012.668230 _12012.483491 | 19613.918138 _19613.888783 | 24930.699029 _24930.689877 | 24936.727305 _ 24936.727305 
CPU time 5s 42s qs 22s 


Example 2. Consider the following sixth-order FSLP: 


D*[y(x)] + $0 a; (a)y (x) + ry(z) = 0, x € (0,5), 


subject to the boundary conditions 
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Figure 1: f-curve according to a = 3.7 


Figure 2: dotted line: Eigenfunction corresponding to 1; dot dashed line: Eigenfunc- 
tion corresponding to 2; dashed line: Eigenfunction corresponding to 3; thick lines: 
Eigenfunction corresponding to 4 of Example 1 


= y'(0) = y (0) = 
=i by =y GO) = 
where 5 < a < 6, q;(x), 0 <j < 5are given by go(z) = —r3(x), q(x) = r(x ); 


qo(x) = r2(x) — rf (x), 93(%) = —2r) (2), qa(z) = —ri(z), r1(z) = 0.0827, 
r(x) = 0.0003x4 — 0.08, oe ) = 10-®x® —0.0014x?. Figure 4 shows the first 
four eigenfunctions of this example. 


Example 3. Consider the following tenth order FSLP: 


—D%y(x) = ry(x) 
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Table 2: Comparison results of first four eigenvalues of Example 2 

Xk a=5.7 a=58 a=59 a= 5.99 
HAM [il] HAM fi HAM fi] HAM [il] 

M1 0.147973 0.144551 0.117816 0.114673 0.094189 0.091268 0.076819 0.074124 

d2 4.575315 4.572390 4.563104 4.560000 4.515779 4.512994 4.461465 4.458254 

X3 44.714848  44.712872 | 45.154858 45.152024 | 46.234038 46.230890 | 47.596156  47.592932 

Na 208.439520 0.699270 | 225.477660 1.826100 | 242.79754 1.923205 | 359.596203 2.595920 
CPU time 304s 308s 327s 307s 


20 fF 


Figure 3: h-curve according to a = 5.99 


Figure 4: dotted line: Eigenfunction corresponding to 1; dot dashed line: Eigenfunc- 
tion corresponding to A2; dashed line: Eigenfunction corresponding to 3; thick lines: 
Eigenfunction corresponding to A4 of Example 2 


subject to the boundary conditions 
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y(0) = (0) = y (0) = y (0) = y™ (0) =0, 
y(n) = 9" (r) = yO (x) = yO (x) = y (a) = 0, 


where 9 < a < 10. The kth exact eigenvalue of this problem is known to be 
0, 


Table 3: Comparison results of first six eigenvalues of Example 3 


Nk a=9.5 a= 9.7 a=99 a= 9.99 a = 9.9999 exact A = k! for a = 10 
mM 1.748599 1.306696 1.075690 1.006707 1.000066 1 
r2 869.080768 974.250744 1018.989895 1023.949824 1024 
x3 46348445365 54048.023602 58515.254210 59043.62485, 59049 
Ya 712914.190485 923496.5577796 1035311.328463 1048443.414051 1048576 
Xs 6378765. 105038 8421530.523672 9620208.712720 9765096.934552 9765625 
Xo 21875877.180515 | 36191783.299387 | 52015671.224970 | 59033451.785699 | 60657173.019026 60466176 
CPU time 73s 78s 77s 70s 81s 31s 


In our proposed method, we use fewer polynomial terms comparing with 
another procedure (see [{11]) and get more accurate results. For instance, 
we just calculate N = 13 and N = 7 terms of approximate solutions, in 
Examples 1 and 2, respectively; While the author in [11] considered N = 18 
and N = 20 terms in aforementioned examples. Moreover, we obtain as the 
same number eigenvalues as Hajji, Al-Mdallal, and Allan reported in [11]. 


5 Conclusion 


In this paper, the HAM has been applied to high order fractional Sturm— 
Liouville problems. The main core of this paper is about to find the eigen- 
values and eigenfunctions of these kind of problems. The proposed method 
has more convenient than the existing methods, since, in spite of using fewer 
terms to approximate the solution, we obtained more accurate numerical re- 
sults. Moreover, it is readily used to solve complex fractional Sturm—Liouville 
problems of an arbitrary high-order. This leads to fewer computations and 
also more efficiency of HAM in comparison with the one proposed in [11]. 
Besides, the HAM on high order FSLPs is validated by verifying its conver- 
gence. 
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